C      ** MAIN PROGRAM FOR HYTRAN PROGRAM  **
C
C
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
      OPEN(5,FILE='HYTRAN.DAT',STATUS='OLD')
      OPEN(6,FILE='HYTRAN.OPT',STATUS='NEW')
C
1     CALL INPUT
      CALL INIT
      GO TO (2,10),IRUN
C
2     CONTINUE
      CALL TMSTEP
      CALL STATE
      CALL BOUDRY
      CALL PHASE
      IF (ITANK1 .NE. 0) CALL TMFUNC
      CALL OUTPUT
      IF (T .LT. TMAX) GO TO 2
      IF (LTAPE .EQ. 0) GO TO 4
      WRITE (6,5) KNTAPE
5     FORMAT (/1X,'KNTAPE=',I10)
4     IF (LPLT .EQ. 0) GO TO 3
      CALL PLOTT
3     CONTINUE
      IF (ISTOP .EQ. 0) GO TO 1
      GO TO 999
C
10    CALL PLOTT
      IF (ISTOP .EQ. 0) GO TO 1
C
999   CLOSE(5,STATUS='KEEP')
      CLOSE(6,STATUS='KEEP')
      STOP
       END
      SUBROUTINE INPUT
C
C     **  ROUTINE TO READ IN DATA  **
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
      NAMELIST/DIN/NJCT,NPIPE,DELT,TMAX,NP,NSORCE,NSINK,NSTANK,HT,C1P,
     1 C2P,C3P,C4P,VOL0,H0,GAMMA,RPC,RPD,VC,H2MIN,VOMAX,D,XL,V0,A,ALPHA,
     2 JJ,IC,C44,ELEV,N,F,NTAU,TIME,TAUV
C
C     PUT PIPE ANGLES IN AS + FOR UP AND - FOR DOWN...COMPUTER WILL CHAN
      READ (5,220) TITLE
220   FORMAT (20A4)
      WRITE (6,3) TITLE
3     FORMAT (1H1,50X,'HYTRAN INPUT DATA'//40X,20A4,//)
      READ(5,DIN)
      WRITE (6,DIN)
C
      READ (5,1) ISTOP,RHO,VARPR,IRUN,ISTART,IDATPE,ITANK1
1     FORMAT (I10,2F10.2,4I10)
      WRITE (6,2)ISTOP,RHO,VARPR,IRUN,ISTART,IDATPE,ITANK1
2     FORMAT (//,2X,'ISTOP',5X,'RHO,LB/FT**3',
     1 8X,'VARPR,FT',8X,'IRUN',4X,'ISTART',4X,'IDATPE',4X,'ITANK1'/
     2 4X,I1,4X,E15.5,5X,E15.5,5X,I1,8X,I1,6X,I5,8X,I1)
C
      IF (ITANK1 .EQ. 0) GO TO 8
C
C     ITANK1=1 FOR SPHERICAL TANK, =2 FOR VETICAL CYLINDRICAL TANK,
C     =3 FOR HORIZONTAL CYLINDRICAL TANK
C
      READ (5,9) RADTK1,XHIGH1,XHIGH2,TKLENG,VFLOIN
9     FORMAT (5F10.3)
      WRITE (6,7) RADTK1,XHIGH1,XHIGH2,TKLENG,VFLOIN
7     FORMAT(/6X,'RADTK1,FT',5X,'XHIGH1,FT',6X,'XHIGH2,FT',6X,
     1 'TKLENG,FT',6X,'VFLOIN,FT**3/S'/4E15.5,E20.5)
8     CONTINUE
C
      READ (5,240)LPRT,ITPRNT,LPLT,ITPLT,LTAPE,ITTPE,TIME2,ITPR2,TIME3,
     1 ITPR3,TIME4,ITPR4
240   FORMAT (6I5,F10.2,I5,F10.2,I5,F10.2,I5)
      WRITE (6,360) LPRT,ITPRNT,LPLT,ITPLT,LTAPE,ITTPE,TIME2,ITPR2,
     1 TIME3,ITPR3,TIME4,ITPR4
360   FORMAT (/1X,'OUTPUT PARAMETERS--LPRT,ITPRNT,LPLT,ITPLT,LTAPE,ITTPE
     1 ,TIME2,ITPR2,TIME3,ITPR3,TIME4,ITPR4'/1X,6I5,F14.2,I5,F14.2,I5,
     2 F14.2,I5)
      IF (LPLT .EQ. 0) GO TO 100
      READ (5,230)NJPLT
230   FORMAT (16I5)
      WRITE (6,380) NJPLT
380   FORMAT (1X,'NJPLT=',I5)
      READ (5,230) (JPLT(I),I=1,NJPLT)
      WRITE (6,390) (JPLT(I),I=1,NJPLT)
390   FORMAT (/1X,'PLOT PARAMETERS--JPLT(I),I=1,2,..'/1X,16I5)
100   CONTINUE
      RETURN
       END
      SUBROUTINE INIT
C
C     ** ROUTINE TO INITIALIZE TRANSIENT PROBLEM **
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
      G=32.2
      ITCTR=0
      KNTAPE=0
      T=0.
      PI=3.14159
C
      IF (NJCT .NE. 0) GO TO 3
51    WRITE (6,4)
4     FORMAT (1X,'SOLUTION TERMINATED DUE TO NJCT=0')
      STOP
3     CONTINUE
      DO114J=1,NPIPE
114   ALPHA(J)=-ALPHA(J)*(3.14159/180.)
      DO2L=1,NJCT
      NL=NP(L)
      DO2M=1,NL
      IF(IC(L,M).EQ.0)GOTO51
      J=JJ(L,M)
      IF(IC(L,M).GT.0)GOTO2
      C2(J)=G/A(J)
      IF (N(J) .EQ. 0) GO TO 51
      XN(J)=N(J)
      DELX(J)=XL(J)/XN(J)
      NN=N(J)+1
      DO53I=1,NN
      V(J,I)=V0(J)
      X=I-1
      H(J,I)=HT(L)-X*DELX(J)*((F(J)*V(J,I)**2)/(2.*D(J)*G))
      VP(J,I)=V(J,I)
      ELJI=ELEV(L)-X*DELX(J)*SIN(ALPHA(J))
53    HP(J,I)=H(J,I)-ELJI
2     CONTINUE
C
      IF (ISTART .EQ. 0) GO TO 5
      READ (ISTART) IDUMP
      DO 6 I=1,IDATPE
      READ (ISTART) T,(HT(L),L=1,NJCT)
      DO 7 J=1,NPIPE
      NN=N(J)+1
7     READ (ISTART) (HP(J,M),VP(J,M),M=1,NN)
6     CONTINUE
5     CONTINUE
C
C     FOR TANK AT JUNCTION 1
      IF (ITANK1 .EQ. 0) GO TO 1
      XDIST=XHIGH1
      IF (ITANK1 .NE. 1) GO TO 9
      VOLTKT=4.*PI*RADTK1**3/3.
      VOLTKG=PI*XHIGH1**2*(RADTK1-XHIGH1/3.)
      GO TO 10
9     IF (ITANK1 .NE. 2) GO TO 11
      VOLTKT=XHIGH2*PI*RADTK1**2
      VOLTKG=XHIGH1*PI*RADTK1**2
      GO TO 10
11    CONTINUE
      THTAT=2.*ATAN(SQRT(2.*XHIGH2*RADTK1-XHIGH2**2)/(RADTK1-XHIGH2))
      VOLTKT=(TKLENG*RADTK1**2/2.)*(THTAT-SIN(THTAT))
      THETAT=2.*ATAN(SQRT(2.*XHIGH1*RADTK1-XHIGH1**2)/(RADTK1-XHIGH1))
      VOLTKG=(TKLENG*RADTK1**2/2.)*(THETAT-SIN(THETAT))
10    CONTINUE
      VOLTKL=VOLTKT-VOLTKG
1     CONTINUE
C
      DO109L=1,NJCT
      HT(L)=HT(L)-ELEV(L)
      IF(NSORCE(L).NE.2)GOTO109
C     ** FOR PUMP HEAD **
      DHEAD(L)=HT(L)
      J1=JJ(L,1)
      N1=N(J1)+1
      SHEAD(L)=H(J1,N1)-ELEV(L)
      HT1(L)=H(J1,N1)
109   CONTINUE
C
      IF (LTAPE .EQ. 0) GO TO 210
      MPTS=(TMAX-T)/(DELT*FLOAT(ITTPE))
      MPTS=MPTS-1
      WRITE (LTAPE) MPTS
210   CONTINUE
C
      CALL OUTPUT
      RETURN
       END
      SUBROUTINE STATE
C
C     ** ROUTINE FOR TRANSIENT CALCULATION **
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
      T=T+DELT
      DO10J=1,NPIPE
      THETA(J)=DELT/DELX(J)
      NN=N(J)
      DO10I=1,NN
      VR(J,I)=(V(J,I+1)-THETA(J)*A(J)*(V(J,I+1)-V(J,I)))/(1.+THETA(J)*(V
     1(J,I+1)-V(J,I)))
      VS(J,I+1)=(V(J,I)-THETA(J)*A(J)*(V(J,I)-V(J,I+1)))/(1.-THETA(J)*(V
     1(J,I)-V(J,I+1)))
      HR(J,I)=H(J,I+1)-THETA(J)*(VR(J,I)+A(J))*(H(J,I+1)-H(J,I))
      HS(J,I+1)=H(J,I)+THETA(J)*(VS(J,I+1)-A(J))*(H(J,I)-H(J,I+1))
10    CONTINUE
C
C     ** FOR INTERIOR POINTS **
      DO11J=1,NPIPE
      NN=N(J)-1
      IF(NN.EQ.0)GOTO11
      DO110I=1,NN
      VP(J,I+1)=0.5*(VR(J,I)+VS(J,I+2)+C2(J)*(HR(J,I)-HS(J,I+2)-DELT*SIN
     1(ALPHA(J))*(VR(J,I)-VS(J,I+2)))-(F(J)*DELT/(2.*D(J)))*(VR(J,I)*ABS
     2(VR(J,I))+VS(J,I+2)*ABS(VS(J,I+2))))
      HP(J,I+1)=0.5*(HR(J,I)+HS(J,I+2)+(A(J)/G)*(VR(J,I)-VS(J,I+2))-DELT
     1*SIN(ALPHA(J))*(VR(J,I)+VS(J,I+2))-(A(J)/G)*(F(J)*DELT/(2.*D(J)))*
     2(VR(J,I)*ABS(VR(J,I))-VS(J,I+2)*ABS(VS(J,I+2))))
110   CONTINUE
11    CONTINUE
C
C     ** FOR CALCULATION OF BOUNDARY POINTS **
      DO12J=1,NPIPE
      I=1
      C1(J)=VS(J,I+1)*(1.+C2(J)*DELT*SIN(ALPHA(J))-F(J)*DELT*ABS(VS(J,I+
     11))/(2.*D(J)))-C2(J)*HS(J,I+1)
      I=N(J)
12    C3(J)=VR(J,I)*(1.-C2(J)*DELT*SIN(ALPHA(J))-F(J)*DELT*ABS(VR(J,I))/
     1(2.*D(J)))+C2(J)*HR(J,I)
C
      RETURN
       END
      SUBROUTINE BOUDRY
C
C     ** ROUTINE FOR BOUNDARY CONDITIONS **
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
      DO13L=1,NJCT
      IF(NSINK(L).NE.0)GOTO27
      IF(IABS(NSORCE(L)) .EQ. 1 .OR. NSORCE(L) .EQ. 3) GO TO 101
      IF(NSORCE(L).EQ.2)GOTO19
C
C     ** FOR NSORCE(L)=0,4,5,ECT. **
      NL=NP(L)
      SUMD2=0.
      SUMDEN=0.
      SUMNUM=0.
      DO14M=1,NL
      J=JJ(L,M)
      XIC(L,M)=IC(L,M)
      IF(IC(L,M).GT.0)GOTO15
      C(L,M)=C1(J)*XIC(L,M)
      GOTO16
15    C(L,M)=C3(J)*XIC(L,M)
16    SUMD2=SUMD2+D(J)**2
      SUMDEN=SUMDEN+C2(J)*D(J)**2
14    SUMNUM=SUMNUM+C(L,M)*D(J)**2
      HT(L)=SUMNUM/SUMDEN
      IF (HT(L) .LE. VARPR) HT(L)=VARPR
      IF(NSTANK(L).NE.0)GOTO23
108   DELQ(L)=0.
26    DO17M=1,NL
      J=JJ(L,M)
      IF(IC(L,M).LT.0)GOTO18
      I=N(J)+1
      HP(J,I)=HT(L)
      VP(J,I)=C3(J)-C2(J)*HP(J,I)-(4.*DELQ(L))/(PI*SUMD2)
      GOTO17
18    I=1
      HP(J,I)=HT(L)
      VP(J,I)=C1(J)+C2(J)*HP(J,I)+(4.*DELQ(L))/(PI*SUMD2)
17    CONTINUE
      GOTO13
C
C     ** SURGE TANK AT JUNCTION   **
23    J=JJ(L,1)
      J1=JJ(L,2)
      HO=H0(L)
      AK2=((PI/4.)*RPD(L)**2*RPC(L))**2
      IF (IC(L,2)) 408,409,410
408   C11=C1(J1)
      C33=C3(J)
      C22=C2(J1)
      C21=C2(J)
      A1=(PI/4.)*D(J)**2
      A2=(PI/4.)*D(J1)**2
      N1=N(J)+1
      GO TO 403
409   C11=-C2(J)*HT(L)
      A2=(PI/4.)*D(J)**2
      C22=C2(J)
      C33=C3(J)
      A1=(PI/4.)*D(J)**2
      N1=N(J)+1
      C21=C2(J)
      GO TO 403
410   C11=C3(J)
      C33=C3(J1)
      C22=C2(J)
      C21=C2(J1)
      A1=(PI/4.)*D(J1)**2
      A2=(PI/4.)*D(J)**2
      N2=N(J1)+1
403   AA=((A1*C21)**2+2.*A1*A2*C21*C22+(A2*C22)**2)/AK2
      BB=2.*(-A1**2*C21*C33-A1*A2*(C22*C33-C21*C11)+A2**2*C11*C22)/AK2-
     1 1.
      CC=((A1*C33)**2-2.*A1*A2*C11*C33+(A2*C11)**2)/AK2+(HO+VARPR)
      IF (HO .LE. (HT(L)-VARPR)) GO TO 404
      BB=BB+2.
      CC=CC-2.*(HO+VARPR)
404   ROOT=SQRT(ABS(BB**2-4.*AA*CC))
      IF (HO .LE. (HT(L)-VARPR)) ROOT=-ROOT
      HT(L)=(-BB+ROOT)/(2.*AA)
      DELH=(HO+VARPR)-HT(L)
      IF (IC(L,2)) 411,416,413
411   VP(J1,1)=C1(J1)+C2(J1)*HT(L)
      VP(J,N1)=C3(J)-C2(J)*HT(L)
      HP(J1,1)=HT(L)
      HP(J,N1)=HT(L)
      GO TO 414
416   VP(J,N1)=C3(J)-C2(J)*HT(L)
      HP(J,N1)=HT(L)
      GO TO 414
413   VP(J,1)=C1(J)+C2(J)*HT(L)
      VP(J1,N2)=C3(J1)-C2(J1)*HT(L)
      HP(J,1)=HT(L)
      HP(J,N2)=HT(L)
414   DELV=(PI/4.)*RPD(L)**2*RPC(L)*SQRT(ABS(DELH))*DELT
      IF (DELH .LE. 0.) DELV=-DELV
      VO=VOL0(L)
      VOL0(L)=VOL0(L)+DELV
      IF (VOL0(L) .LE. 0.) VOL0(L)=0.00001
      H0(L)=H0(L)*(VO/VOL0(L))**GAMMA(L)
      IF (H0(L) .GT. H2MIN(L)) GO TO 415
      H0(L)=H2MIN(L)
      VOL0(L)=VOMAX(L)
415   GO TO 13
C
C     PUMP BETWEEN TWO PIPES OF EQUAL DIAMETER
19    J=JJ(L,1)
      J1=JJ(L,2)
      I=N(J)+1
      IABORT=19
      IF((C1P(L).EQ.0.).AND.(C3P(L).EQ.0.))GOTO51
      IF((C2P(L).EQ.0.).AND.(C3P(L).EQ.0.))GOTO51
      IF(C3P(L).EQ.0.)GOTO51
      AA=C1P(L)
      BB=C2P(L)-1./C2(J)-1./C2(J1)
      CC=C3(J)/C2(J)+C1(J1)/C2(J1)+C3P(L)
      VP(J,I)=(-BB-SQRT(ABS(BB**2-4.*AA*CC)))/(2.*AA)
      IF(VP(J,I).LE.0.)VP(J,I)=0.
      VP(J1,1)=VP(J,I)
      HP(J,I)=(VP(J,I)-C3(J))/(-C2(J))
      HP(J1,1)=(VP(J1,1)-C1(J1))/C2(J1)
      SHEAD(L)=HP(J,I)-ELEV(L)
      DHEAD(L)=HP(J1,1)-ELEV(L)
      HT(L)=HP(J1,1)
      GOTO13
C
C     PUMP OR RESERVOIR AT END OF PIPE
101   M=1
      J=JJ(L,M)
      I=1
      IF (NSORCE(L) .EQ.(-1)) I=N(J)+1
      IABORT=101
      IF((C1P(L).EQ.0.).AND.(C3P(L).EQ.0.))GOTO51
      IF((C2P(L).EQ.0.).AND.(C3P(L).EQ.0.))GOTO51
      IF(C3P(L).EQ.0.)GOTO51
      IF(C1P(L).EQ.0.)GOTO20
      AA=C1P(L)*C2(J)
      BB=C2P(L)*C2(J)-1.
      CC=C3P(L)*C2(J)+C1(J)
      VP(J,I)=(-BB-SQRT(BB**2-4.*AA*CC))/(2.*AA)
      IF(VP(J,I).LE.0.)VP(J,I)=0.
      GOTO22
20    IF(C2P(L).EQ.0.)GOTO21
      VP(J,I)=(C3P(L)-(C1(J)/C2(J)))/(C2P(L)+1./C2(J))
      IF(VP(J,I).LE.0.)VP(J,I)=0.
      GOTO22
21    CONTINUE
      VP(J,I)=C1(J)+C2(J)*C3P(L)
      IF(NSORCE(L) .EQ. (-1))VP(J,I)=C3(J)-C2(J)*C3P(L)
      IF(NSORCE(L).EQ.3)GOTO111
C
C     FOR PUMP BETWEEN TWO PIPES THIS IS HEAD INCREASE
22    HT(L)=C1P(L)*VP(J,I)**2+C2P(L)*VP(J,I)+C3P(L)
      IF (HT(1) .LT. 0.) HT(1)=0.
      HP(J,I)=HT(L)
      GOTO13
C
C     RESERVOIR WITH CHECK VALVE TO PREVENT BACK FLOW ... ALSO SIMULATES
C     HYDROPNEUMATIC TANK WITH CONSTANT GAS BLANKET PRESSURE
111   IF(VP(J,I).LE.0.)VP(J,I)=0.
      HP(J,I)=(VP(J,I)-C1(J))/C2(J)
      HT(L)=HP(J,I)
      GOTO 13
C
C     ** FOR VALVES  **
27    NV=NSINK(L)
      M=1
      J=JJ(L,M)
      I=N(J)+1
      IPASS=1
      IABORT=27
      IF(NV .GT. 2) GO TO 81
      NTAB=NTAU(L)
      DO 80 K=1,NTAB
      TMD(K)=TIME(K,L)
      TAUD(K)=TAUV(K,L)
80    CONTINUE
      CALL INTP(T,TAU,NTAB,TMD,TAUD)
      GO TO 28
81    IF (NV .GT. 7) GO TO 51
      TAU=1.
28    C4=C44(L)  *TAU**2
      GO TO    (37,36,37,36,37,37,36),NV
C
C     VALVE OR ORIFICE AT END OF PIPE DISCHARGING TO ATMOSPHERE
36    VP(J,I)=-C4/2.+SQRT(ABS((C4/2.)**2+C3(J)*C4))
      HP(J,I)=(C3(J)-VP(J,I))/C2(J)
      IF (HP(J,I) .LE. ELEV(L)) GO TO 116
      IF (NV .EQ. 7) GO TO 104
      GO TO 105
116   HP(J,I)=ELEV(L)
      VP(J,I)=C3(J)-C2(J)*ELEV(L)
105   HT(L)=HP(J,I)
      GOTO13
C
C     PRESSURE RELIEF VALVE
104   XVALVE=C1P(L)
      IF(HP(J,I).LE.XVALVE)VP(J,I)=0.
      HP(J,I)=(C3(J)-VP(J,I))/C2(J)
      GOTO105
C
C     VALVE OR ORIFICE BETWEEN TWO PIPES OF EQUAL DIAMETER
37    J1=JJ(L,M+1)
      I1=1
      C1PC3=C1(J1)+C3(J)
      IF((C1PC3.GE.0.).AND.(TAU.GT..00001))GO TO 38
      IF (NV .EQ. 5) TAU=0.0
      IF(TAU.LT..00001)GO TO 39
C
C     NEGATIVE FLOW THRU VALVE OR ORIFICE
      VP(J,I)=C4-SQRT(C4**2-C4*C1PC3)
48    VP(J1,I1)=VP(J,I)
      HP(J,I)=(C3(J)-VP(J,I))/C2(J)
      HP(J1,I1)=(VP(J1,I1)-C1(J1))/C2(J1)
      HT(L)=HP(J,I)
      HT1(L)=HP(J1,I1)
      IF (IPASS .EQ. 2) GO TO 13
C
C     C1P(L) BELOW IS HEAD AT WHICH PRESSURE REDUCING VALVE IS ACTIVATED
      IF((NV .EQ. 6) .AND. (HP(J1,I1) .GE. C1P(L))) GO TO 115
      GOTO13
C
C     POSITIVE FLOW THRU VALVE OR ORIFICE
38    VP(J,I)=-C4+SQRT(C4**2+C4*C1PC3)
      GOTO48
39    VP(J,I)=0.
      GO TO 48
115   IPASS=2
C
C     C2P(L) BELOW IS HEAD AT WHICH PRESSURE REDUCING VALVE IS COMPLETEL
C     CLOSED
      TAU=((C2P(L)-HP(J1,I1))/(C2P(L)-C1P(L)))**2
      IF (HP(J1,I1) .GE. C2P(L)) TAU=0.
      C4=C44(L)  *TAU**2
      GO TO 37
13    CONTINUE
      RETURN
C
51    WRITE (6,52) IABORT
52    FORMAT (1X,'SOLUTION TERMINATED DUE TO BAD DATA'/
     1 1X,'IABORT=',I5,2X,'IN BOUNDARY')
      STOP
      END
      SUBROUTINE PHASE
C
C     ** ROUTINE FOR PHASE **
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
      DO40L=1,NJCT
      NL=NP(L)
      DO112M=1,NL
      IF(IC(L,M).GT.0)GOTO112
      J=JJ(L,M)
      NN=N(J)+1
      DO 115 I=1,NN
      X=I-1
      ELJI=ELEV(L)-X*DELX(J)*SIN(ALPHA(J))
      VAPHJI=ELJI+VARPR
      IF(HP(J,I).LE.VAPHJI)GOTO500
      GOTO501
C
C      LIQUID FLOW WITH GAS BUBBLE
500    CONTINUE
      XY(J,I)=1.
      HP(J,I)=VAPHJI
      IF(I.EQ.NN)GOTO501
      IF (I .EQ. 1) GO TO 501
      VCP=VR(J,I)*(1.-C2(J)*DELT*SIN(ALPHA(J))-F(J)*DELT*ABS(VR(J,I))/(2
     1.*D(J)))+C2(J)*HR(J,I)
      VCM=VS(J,I+1)*(1.+C2(J)*DELT*SIN(ALPHA(J))-F(J)*DELT*ABS(VS(J,I+1)
     1)/(2.*D(J)))-C2(J)*HS(J,I+1)
      VP(J,I)=VCP-C2(J)*HP(J,I)
      VPP(J,I)=VCM-C2(J)*HP(J,I)
      XCAV(J,I)=XCAV(J,I)+(VPP(J,I)-VP(J,I))*DELT
      IF(XCAV(J,I).LE.0.)GOTO502
      GOTO501
502   XY(J,I)=0.
      VP(J,I)=0.5*(VCM+VCP)
      HP(J,I)=VAPHJI+(A(J)/(2.*G))*(VCP-VCM)
      XCAV(J,I)=0.
C
501   V(J,I)=VP(J,I)
      H(J,I)=HP(J,I)
      HP(J,I)=HP(J,I)-ELJI
      IF (HP(J,I) .LE. VARPR) HP(J,I)=VARPR
115   CONTINUE
112   CONTINUE
40    CONTINUE
      DO50L=1,NJCT
      VAPHTL=ELEV(L)+VARPR
      HT(L)=HT(L)-ELEV(L)
      IF (HT(L) .LE. VAPHTL) HT(L)=VARPR
50    CONTINUE
      RETURN
      END
      SUBROUTINE OUTPUT
C
C     ** ROUTINE TO OUTPUT TRANSIENT SOLUTION **
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
      IF (LPRT .EQ. 0) GO TO 10
      IF (TIME2 .LE. 0.0) GO TO 77
      IF (T .GT. TIME2) ITPRNT=ITPR2
17    IF (TIME3 .LE. 0.0) GO TO 77
      IF (T .GT. TIME3) ITPRNT=ITPR3
      IF (TIME4 .LE. 0.0) GO TO 77
      IF (T .GT. TIME4) ITPRNT=ITPR4
77    CONTINUE
      IF (MOD(ITCTR,ITPRNT) .NE. 0) GO TO 10
      WRITE (6,4)T
4     FORMAT(1H1,' SOLUTION OF FLUID SYSTEM TRANSIENT AT TIME=',E15.5,
     1 2X,'SEC',//)
      WRITE (6,5)
5     FORMAT (1X,'PIPE NO.',6X,'PRESSURE HEAD UP STREAM END',4X,
     1 'PRESSURE HEAD DN STREAM END',10X,'VEL. UP STREAM END',4X,
     2 'VEL. DN STREAM END'/20X,'FT.',10X,'PSI',14X,'FT.',11X,'PSI',
     3 23X,'FT/SEC',16X,'FT/SEC',/)
      DO 8 J=1,NPIPE
      I=1
      NN=N(J)+1
      HPIPSI=HP(J,I)*RHO/144.
      HPNPSI=HP(J,NN)*RHO/144.
8     WRITE (6,9) J,HP(J,I),HPIPSI,HP(J,NN),HPNPSI,VP(J,I),VP(J,NN)
9     FORMAT (3X,I4,5X,4E15.5,10X,E15.5,5X,E15.5)
C
      IF (ITANK1 .EQ. 0) GO TO 1
      IF (T .LE. 1.0E-10) GO TO 18
      WRITE (6,19)
19    FORMAT (/1X,'DATA AT NJCT=1 (TANK) HAVE BEEN ADJUSTED',/)
18    WRITE (6,2) VOLTKT,VOLTKG,VOLTKL,C3P(1),ELEV(1),XL(1),DELX(1)
2     FORMAT (/1X,'VOLTKT,VOLTKG,VOLTKL,C3P(1),ELEV(1),XL(1),DELX(1) ARE
     1,RESPECTIVELY'/10X,7E15.5)
1     CONTINUE
C
      DO102L=1,NJCT
102   IF(NSORCE(L).EQ.2)WRITE(6,103)L,SHEAD(L),DHEAD(L)
103   FORMAT(1H0,16HPUMP AT JUNCTION,I3,5X,13HSUCTION HEAD=,F8.1,5X,15HD
     1ISCHARGE HEAD=,F8.1)
C
10    IF (LTAPE .EQ. 0) GO TO 40
      IF (MOD(ITCTR,ITTPE) .NE. 0) GO TO 40
      WRITE (LTAPE) T,(HT(L),L=1,NJCT)
      DO 7 J=1,NPIPE
      NN=N(J)+1
7     WRITE (LTAPE) (HP(J,I),VP(J,I),I=1,NN)
      KNTAPE=KNTAPE+1
C
40    CONTINUE
      ITCTR=ITCTR+1
      RETURN
       END
      SUBROUTINE TMSTEP
C
C     ** ROUTINE TO CHECK TIME STEP FOR STABILITY **
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
      TOL=1.0E-6
      DTMIN=DELX(1)/(V(1,1)+A(1))
      DO 1 J=1,NPIPE
      NN=N(J)+1
      DO 2 I=1,NN
      DTEST=DELX(J)/(V(J,I)+A(J))
      IF (DTMIN .GT. DTEST) DTMIN=DTEST
      IF (DTMIN .GT. TOL) GO TO 2
      WRITE (6,9) J,I,DELX(J),V(J ,I),A(J),DTMIN
9     FORMAT (/1X,'TIME STEP TOO SMALL---J,I,DELX(J),V(J,I),A(J),DTMIN'/
     11X,2I10,4E20.5)
      LPRT=1
      ITPRNT=1
      CALL OUTPUT
      STOP
2     CONTINUE
1     CONTINUE
8     IF (DELT .GT. DTMIN) GO TO 4
      RETURN
4     DELT=DELT/2.
      ITPRNT=2*ITPRNT
      ITPLT=2*ITPLT
      ITTPE=2*ITTPE
      ITPR2=2*ITPR2
      ITPR3=2*ITPR3
      ITPR4=2*ITPR4
      WRITE (6,5) DELT
5     FORMAT (/1X,'TIME STEP HAS BEEN HALVED,DELT=',E15.6)
      GO TO 8
      END
      SUBROUTINE  TMFUNC
C
C     ROUTINE FOR TANK AT JUNCTION ONE
C
      DIMENSION XC(11),FX(11),XI(18),FD(18)
C
      COMMON /IN1/ NP(30),NSORCE(30),NSINK(30),NSTANK(30),HT(30),
     1 C1P(30),C2P(30),C3P(30),VOL0(30),H0(30),GAMMA(30),RPC(30),
     2 RPD(30),VC(30),H2MIN(30),VOMAX(30),C4P(30)
      COMMON /INA/ D(40),XL(40),V0(40),A(40),ALPHA(40),F(40)
      COMMON /IN2/C44(30),IC(30,5),JJ(30,5),NTAU(30),TIME(10,30),
     1 TAUV(10,30),TMD(10),TAUD(10)
      COMMON /IN3/ NJCT,NPIPE,DELT,T,TMAX,ITPRNT,ISTOP,LPRT,LPLT,
     1 ITPLT,LTAPE,ITTPE,NJPLT,JPLT(30),TIME2,ITPR2,TIME3,ITPR3
     2 ,TIME4,ITPR4
      COMMON /BK1/ XN(40),DELX(40),THETA(40),N(40),C1(40),C3(40),C2(40),
     1 DELQ(30),HT1(30),DHEAD(30),SHEAD(30),ELEV(30)
      COMMON /BK2/ V(40,50),H(40,50),VS(40,50),HS(40,50),VR(40,50),
     1 HR(40,50),VP(40,50),HP(40,50),XIC(30,5),C(30,5)
      COMMON /BK3/ G,PI,RHO,VARPR,TITLE(20),ITCTR,KNTAPE,IRUN
     1 ,ISTART,IDATPE
      COMMON /BUBBL/ XY(40,5),VPP(40,50),XCAV(40,50)
      COMMON /TANK/ ITANK1,RADTK1,XHIGH1,XHIGH2,XDIST,VOLTKT,VOLTKG,
     1 VOLTKL,VFLOIN,TKLENG
C
C     DATA FOR SPHERICAL TANK
      DATA KN,XC,FX/11,0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,0.,
     1   0.028,0.104,0.216,0.352,0.5,0.648,0.784,0.896,0.972,1.0/
C
C     DATA FOR HORIZONTAL CYLINDRICAL TANK
      DATA KM,XI,FD /18,0.,0.01,0.08,0.27,0.6,1.09,1.72,2.47,3.14,3.26,
     1 4.04,4.76,5.35,5.8,6.08,6.23,6.26,6.28, 0.,0.4,0.8,1.2,1.6,2.0,
     2 2.4,2.8,3.1416,3.2,3.6,4.0,4.4,4.8,5.2,5.6,6.0,6.28/
C
C     ITANK1=1 FOR SPHERICAL TANK, =2 FOR VETICAL CYLINDRICAL TANK,
C     =3 FOR HORIZONTAL CYLINDRICAL TANK
C
      FLOUT=VP(1,1)*DELT*PI*D(1)**2/4.-VFLOIN
      IF (XDIST-XHIGH2)4,3,3
4     VOLTKG=VOLTKG+FLOUT
      VOLTKL=VOLTKT-VOLTKG
      IF (ITANK1 .NE. 1) GO TO 5
      VRATIO=VOLTKG/VOLTKT
      CALL INTP (VRATIO,XANS,KN,FX,XC)
      XDISTO=XDIST
      XDIST=2.*RADTK1* XANS
      DDIST=XDIST-XDISTO
      GO TO 6
5     IF (ITANK1 .NE. 2) GO TO 7
      DDIST=FLOUT/(PI*RADTK1**2)
      XDIST=XDIST+DDIST
      GO TO 6
7     CONTINUE
      VRATIO=2.*VOLTKG/(TKLENG*RADTK1**2)
      CALL INTP (VRATIO,XANS,KM,XI,FD)
      XDISTO=XDIST
      XDIST=RADTK1*(1.-COS(0.5*XANS))
      DDIST=XDIST-XDISTO
6     CONTINUE
      C3P(1)=C3P(1)-DDIST
      RETURN
C
3     DDIST=VP(1,1)*DELT
      XL(1)=XL(1)-DDIST
      XN(1)=N(1)
      DELX(1)=XL(1)/XN(1)
      ELEV(1)=ELEV(1)-DDIST*SIN(ALPHA(1))
      C3P(1)=C3P(1)-DDIST*SIN(ALPHA(1))
      RETURN
      END
      SUBROUTINE INTP (X,Y,NN,A,B)
      DIMENSION A(1),B(1)
C
C     **  ROUTINE FOR LINEAR INTERPOLATION  **
C     ** ARRAY A MUST BE IN ASCENDING ORDER **
C
C     TO CHECK IF X IS WITHIN TABULATED RANGE
C     SET Y EQUAL TO END-VALUE IF X IS OUT OF RANGE
C
      IF (X .GT. A(1)) GO TO 10
      Y=B(1)
      RETURN
10    IF (X .LT. A(NN)) GO TO 20
      Y=B(NN)
      RETURN
20    CONTINUE
C
C     TO FIND X1 AND Y1
C
      DO 30 I=1,NN
      IM1=I-1
      TEST=X-A(NN-IM1)
      IF (TEST) 30,34,38
38    X1=A(NN-IM1)
      Y1=B(NN-IM1)
      GO TO 40
34    Y=B(NN-IM1)
      RETURN
30    CONTINUE
40    CONTINUE
C
C    TO FIND X2 AND Y2
C
      DO 50 I=1,NN
      TEST=X-A(I)
      IF (TEST .GT. 0.) GO TO 50
      X2=A(I)
      Y2=B(I)
      GO TO 60
50    CONTINUE
60    CONTINUE
C
C    TO PERFORM LINEAR INTERPOLATION
C
      Y=Y1-(Y1-Y2)*(X1-X)/(X1-X2)
      RETURN
      END
      SUBROUTINE PLOTT
C
C     ** ROUTINE TO PLOT TRANSIENT SOLUTIONS ON SC4460 **
C
      RETURN
      END
